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Abstract 

We show that by requiring positivity of the longitudinal pressure it is possible to constrain the 
initial conditions one can use in 2nd-order viscous hydrodynamical simulations of ultrarelativistic 
heavy- ion collisions. We demonstrate this explicitly for 0+1 dimensional viscous hydrodynamics 
and discuss how the constraint extends to higher dimensions. Additionally, we present an analytic 
approximation to the solution of 0+1 dimensional 2nd-order viscous hydrodynamical evolution 
equations appropriate to describe the evolution of matter in an ultrarelativistic heavy-ion collision. 
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I. INTRODUCTION 



The main goal of experiments which perform ultrarelativistic heavy-ion colhsions is to 
produce and study the properties of a deconfined plasma of quarks and gluons. This new 
state of matter, the quark-gluon plasma (QGP), is expected to be formed once the temper- 
ature of nuclear matter exceeds a critical temperature of Tq ~ 200 MeV. Such experiments 
have already been underway for nearly a decade at the Relativistic Heavy Ion Collider 
(RHIC) and higher-energy runs are planned at Large Hadron Collider (LHC). Historically, 
in order to make phenomenological predictions for experimental observables, fluid hydro- 
dynamics has been used to model the space-time evolution and non-equilibrium properties 
of the expanding matter. For the description of nuclear matter by fluid hydrodynamics to 
be valid the microscopic interaction time scale must be much shorter than the macroscopic 
evolution time scale. However, the hot and dense matter created in these experiments is 
rather small in transverse extent and expands very rapidly causing the range of validity of 
hydrodynamics to be limited. 

After the first results of RHIC, it was somewhat of a surprise that ideal hydrodynamics 
could reproduce the hadron transverse momentum spectra in central and semi-peripheral 
collisions. This included their anisotropy in non-central collisions which is measured by 
the elliptic flow coefficient, V2{pt)- Ideal hydrodynamical models were fairly successful in 
describing the dependence of V2 on the hadron rest mass for transverse momenta up to 
about 1.5-2 GeV/c [H, Q, [s], 0] . This observation led to the conclusion that the QGP formed 
at RHIC could have a short thermalization time (tq < 1 fm/c) and a low shear viscosity. 
As a result it was posited that the matter created in the experiment behaves like a nearly 
perfect fluid starting at very early times after the collision. However, recent results from 
viscous hydrodynamical simulations which include all 2nd-order transport coefficients con- 
sistent with conformal symmetry [sl have shown that estimates of the thermalization time 
are rather uncertain due to poor knowledge of the proper initial conditions, details of plasma 
hadronization, subsequent hadronic cascade, etc.^. As a result, it now seems that thermaliza- 
tion times of up to tq ~ 2 fm/c are not completely ruled out by RHIC data. Faced with this 
challenge it has been recently suggested that it may be possible to experimentally constrain 
To by makin g us e of high-energy electromagnetic probes such as dileptons jl, 10, [ll, 13] and 



photons [131, M, [15 



As mentioned above, one of the key ingredients necessary to perform any numerical sim- 
ulation using fluid hydrodynamics is the proper choice of initial conditions at the initially 
simulated time (tq). These initial conditions include the initial fluid energy density e, the 
initial components of the fluid velocity and the initial shear tensor U^'^. Once the set 
of initial conditions is known, it is "simple" to follow the subsequent dynamics of the fluid 
equations in simulations. At the moment there is no first principles calculation that allows 
one to determine the initial conditions necessary. Two different approaches are currently 
used for numerical simulations of fluids in heavy-ion collisions: Glauber type [l6| or Colored- 
Glass-Condensate (CGC) initial conditions^. The uncertainty in the initial conditions in- 
troduces a systematic theoretical uncertainty when, for example, the transport coefficient 



^ For more about the application of viscous hydrodynamics to heavy-ion phenomenology we refer the reader 

to Rcfs. B,flB[3. 

^ For a recent review on the initial conditions based on the CGC approach see Ref. [17| and references 
therein. 
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7]/ s is extracted from experimental data [sl, H, 0, [!]• This is due to the fact that when the 
initial energy density profile is fixed using CGC-based initial conditions (isl . ITol . 20| . one 
obtains larger initial spatial eccentricity and momentum anisotropy when compared with 
the Glauber model. Moreover, the values of the components of the shear tensor 11^^ at tq 
are also affected by the choice of either CGC or Glauber initial conditions (see discusion in 
Sect. 4 of Ref. [2l|). In the case of Glauber initial conditions the shear tensor is completely 
unconstrained. In the case of CGC initial conditions there is a prescription for calculating 
the initial shear; however, with CGC initial conditions the longitudinal pressure is zero due 
to the assumption of exact boost invariance and the subsequent thermalization of the system 
could completely change the initial shear obtained in the CGC approximation. Therefore, 
in both cases it would seem that the initial shear is completely unconstrained. 

Given these uncertainties it would be useful to have a method which can help to constrain 
the allowed initial conditions used in hydrodynamical simulations. In this work we derive 
general criteria which impose bounds on the initial time tq at which one can apply 2nd- 
order viscous hydrodynamical modeling of the matter created in ultrarelativistic heavy-ion 
collision. We do this firstly by requiring the positivity of the effective longitudinal pressure 
and secondly by requiring that the shear tensor be small compared to the isotropic pres- 
sure. Based on these requirements we find that, for a given set of transport coefficients, 
the allowed minimum value of tq is non-trivially related with the initial condition for the 
shear tensor, n'^'^(ro) = IIq'^, and the energy density e(ro) = en. To make this explicit we 



study 0-1-1 dimensional 2nd-order viscous hydrodynamics [22, l23|, |24| . where the transport 



coefficients are either those of a weakly-coupled transport theory j25l . 26. 27| or those ob 



tained from a strongly-coupled A/" = 4 supersymmetric (SYM) plasma 23|, |2J]. We then 
show how the constraints derived from the 0-1-1 dimensional case can be used to estimate 
where higher dimensional simulations will cease to be physical/trustworthy. Our technique 
is complementary to the approach of Molnar and Huovinen [28] which uses kinetic theory to 
assess the applicability of hydrodynamics. In contrast to their work, here we do not invoke 
any other physics other than hydrodynamical evolution itself and merely require that it be 
reasonably self-consistent. 

The work is organized as follows: in Sec. [Ill we review the basic setup of 2nd-order viscous 
hydrodynamics formalism and its application to a -|- 1 dimensional boost invariant QGP 
(either in the weakly or strongly coupled limits). In Sec. IIIII we present an approximate 
analytical solution to the equations of motion for a -|- 1 dimensional system. In Sec. IIV[ 
we present our analytical and numerical results in both the strong and weak coupling limits 
of the 0-1-1 dimensional QGP. In Sec. [Vl we present our conclusions. 

II. BASIC SETUP 

In this section we briefly review the general framework of 2nd-order viscous hydrody- 
namics equations for a conformal fluid, i.e. we will consider just shear viscosity and neglect 
bulk viscosity. We will also ignore heat conduction. The energy-momentum tensor for a 
relativistic fluid in the presence of shear viscosity is given by^: 

T^"" = eu^'u'' -p^^"" + (2.1) 



The notation we use along the text is summarized in the Appendix [X] 
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where e and p are the fluid energy density and pressure, is the normahzed fluid four- 
velocity {u^u^ = 1) and U'^" is the shear tensor which has two important properties: (1) 
= and (2) Uf^U'^'' = 0. Requiring conservation of energy and momentum, D^T^'^ = 0, 
gives the space-time evolution equations for the fluid velocity and the energy density: 

De = -(e + VX + l^^'^iuu,) , (2.2) 

where D^^ is the geometric covariant derivative, D = u'^Da is the comoving time derivative 
in the fluid rest frame and = A'^'^Da is the spatial derivative in the fluid rest frame. 
The brackets ( ) construct terms which are symmetric, traceless, and orthogonal to the fluid 
velocity (see Appendix [Al for its deflnition). 

To obtain a complete solvable system of equations viscous hydrodynamics requires an 
additional equation of motion for the shear tensor. This is accomplished by expanding the 
equations of motion to second order in gradients. It has been found that at zero-chemical 



potential in a conformal fluid in any curved space-time, the shear tensor satisfles [23|, |24| |: 



4, 
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.u 



- An<^n->^ + ^^<^^^>^ - ^uj<'^xio''>' , (2.3) 

where u^i, = — Vj^Mj/J is a symmetric operator that represents the fluid vorticity and W^^^^l^ 
and R^^ are the Riemann and Ricci tensors, respectively. The coefficients t^^, k,\i^\2 and 
A3 are the transport coefficients required by conformal symmetry. 



A. 0-|-l Dimensional Conformal 2nd-Order Viscous Hydrodynamics 



Let us consider a system expanding in a boost invariant manner along the longitudinal 
(beamline) direction with a uniform energy density along the transverse plane. For this 
simplest heavy-ion collision model, it is enough to consider expansion in a flat space. Also 
for this simple model, there is no fluid vorticity, and the energy density, the shear viscous 
tensor and the fluid velocity only depend on proper time r. For this 0+1 dimensional 
model the 2nd-order viscous hydrodynamic equations (Eqs. 02.21) and (12.31) ) are rather simple 

roper time, r = ^/t"^ — , and space-time rapidity, 
23: 



in the conformal limit. In terms of 
C = arctanh(2;/t), these are given by [2 



e -Fp _^ n 

r r 



4?7 



n 



3r 



Al 



2r^ri' 



(2.4) 
(2.5) 



where e is the fluid energy density, p is the fluid pressure, 11 = 11^ is the CC component of 
the fluid shear tensor, ri is the fluid shear viscosity, r,r is the shear relaxation time, and Ai 
is a coefficient which arises in complete 2nd-order viscous hydrodynamical equations either 
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Transport coefficient 


Weakly-coupled QCD 


Strongly-coupled M = 4 SYM 


f] = rj/s 


~ 1/(5'^ logs-) 


1/(4^) 




6r?/r 


(2-Iog2)/(2^r) 


Ai 


(4.1 ^ 5.2) fs/T 


2 ?72s/r 



TABLE I: Typical values of the transport coefficients for a weakly-coupled QGP 2^, 26, 27 1 and a 
strongly coupled M = 4 SYM plasma |23l.l24|. 



in the strong 23, 24] or weakly coupled limit [13, 25, 26, 27, 2^. The Navier-Stokes limit is 
recovered upon taking — » and Ai — >• in which case one obtains IlNavier-stokes = '^vl (3t). 

These coupled differential equations are completed by a specification of the equation 
of state which relates the energy density and the pressure through p = p{e) and initial 
conditions. For 0+1 dimensional dynamics one must specify the energy density and 11 at 
the initial time, eo = e(ro) and Ho = n(ro), where tq is the proper-time at which one begins 
to solve the differential equations. 



B. Specification of equation of state and dimensionless variables 

In the following analysis we will assume an ideal equation of state, in which case we have 

P = , (2.6) 

where for quantum chromodynamics with Nc colors and Nf quark flavors, A'dof = 2(A'"^ — 
1) + 7NcNf/2 which for Nc = 3 and Nf = 2 is N^of = 37. The general method used below, 
however, can easily be extended to a more realistic equation of state. 

In the conformal limit the trace of the four- dimensional stress tensor vanishes requiring 
e = 3p which, using Eq. (12.61) . allows us to write compactly 

30 



6 = (T/7)^ with 7 = • (2.7) 

Likewise we can simplify the expression for the entropy density, s, using the thermodynamic 
relation Ts = e + p to obtain s = 4e/3T or equivalently 

s = ^e'/\ (2.8) 
37 

When solving Eqs. (12. 4p and (12.51) it is important to recognize that the transport coeffi- 
cients depend on the temperature of the plasma and hence on proper-time. We summarize 
in Table [T] the values of the transport coefficients in the strong and weak coupling limits. 
We point out that these are not universal relations as explained below in Sees. Ill Gl and III D[ 
The reader should note that in either the strong or weak coupling limit the coefficients 
and Ai are proportional to oc and Ai oc ffs/T. This suggests that we can parametrize 
both coefficients as: 

r. = Y^ (2-9a) 

^2^ s 



Ai = c,,r/^(-), (2.9b) 
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where we have introduced a dimensionless version of the shear viscosity 



T] = rj/s . 



(2.10) 



In our analysis we assume that fj is independent of time.'^ The dimensionless numbers f], 
and carry all of the information about the particular coupling limit we are considering. 

Using the ideal gas equation of state [Eqs. (12. 7p and (12. Sp ]. the parametrization (12. 9p of 
r,r and Ai can be rewritten in terms of the energy density e: 



7e 
i 

372 



1/4 ' 



Ai = ^ cai r] e 



2^1/2 



(2.11a) 
(2.11b) 



To remove the dimensionful scales and rewrite the fluid equations in a more explicit form 
we define the following dimensionless variables: 



n = n/eo 



(2.12a) 
(2.12b) 
(2.12c) 



where tq is the proper-time at which the hydrodynamic evolution equations start to be 
integrated and eo is the energy density at tq. 

After replacing the dimensionless variables (I2.12p in the parametrization (12.111) and 
Eqs. (12. 4p and (12. 5p . we rewrite the fluid equations: 



Tdfe + - e 

n+ 



7A;eV4 



n = 0, 



3 r 



-3/4 



16 77 e 

97 A; f 



+ 



3CA, 



where k = Toe^^. Note that in terms of (I2.12p the boundary conditions are specified at 
f = 1 where e = 1 and n(f = 1) = IIo which is a free parameter. When the hydrodynamical 
equations are written in the form given above [Eq. (12.131) ] all information about the initial 
proper-time and energy density is encoded in the parameter k and all information about the 
equation of state is encoded in the parameter 7. 



(2.13a) 
(2.13b) 



C. Strong coupling limit 

Motivated and guided by the AdS/CFT correspondence Baier et. al j23[ and the Tata 
group 2J] have recently shown that new transport coefficients arise in a complete theory of 
second order relativistic viscous hydrodynamics. They also estimate their values at infinite 
t'Hooft coupling for A/" = 4 SYM theory at finite temperature. Different calculations for a 



Including a temperature-dependent shear viscosity does not change our observations fundamentally; how- 
ever, there will be quantitative effects which will be elaborated upon in a forthcoming publication. 
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finite t'Hooft coupling within the same theory have been carried out 0, [m, S H S, M 
A remarkable aspect is that, while at first the strong t'Hooft coupling limit of the transport 



coefficients was expected to be universal [36|, l37| , there is now evidence that these coefficients 



are not universal [38|, |39|, |40|, |4l|, |42| . Faced with this complication one is forced to make 
a choice as to which dual theory to consider. Here we will consider the values obtained in 
A/" = 4 SYM at infinite t'Hooft coupling as used in jH, 2^ as our typical strong coupling 
values. One can expect that these coefficients change in strongly-coupled QCD compared to 
A/" = 4 SYM theory at infinite t'Hooft limit. Nevertheless, we take these values over from 
strongly-coupled A/" = 4 SYM in order to get a feeling for what to expect in this regime. 

Expressed in terms of the dimensionless transport coefficients defined above, typical values 
of the strongly coupled transport coefficients are 



^ 47r 



Ctj- 



2 - log 2 
2^ 



(2.14) 



D. Weak coupling limit 

Contrary to the case of A/" = 4 SYM at infinite coupling, in the case of QCD, where 
there is a running coupling and inherent scale dependence, the various transport coefficients 
are not fixed numbers but instead depend on the renormalization scale. In this limit the 



transport coefficients necessary have been calculated completely to leading order [25|, l26|, |27 
Higher order corrections to some transport coefficients from finite-temperature perturbation 
theory show poor convergence 43|, |4^ which is similar to the case for the thermodynamical 
potential; however, resummation techniques can dramatically extend the range of conver- 
gence of finite-temperature perturbation theory in the case of static quantities and can, in 
the future, also be applied to dynamical quantities^. Until such resummation schemes are 
carried out for dynamical quantities, the values of the leading-order weak-coupling trans- 
port coefficients in Table [J can only be considered as rough guides to the values expected 
phenomenologically. Using this rough guide the value of f] from finite-temperature QCD 

^ 1 at realistic couplings ((7 ~ 2 ^ 3). In this work we 
10 / (Air) in the weakly-coupled limit in order to compare 
In our analysis for the weak coupling 



calculations [2n, 27 



is 77/3 ~ 0.5 
will assume a typical value of fj = 
with the results obtained in the strong coupling limit 
limit, we will use 



10 

= 6?7 , 
9 



(2.15) 



See Ref. 



45[ and references therein. 
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E. Momentum space anisotropy 



We introduce the dimensionless parameter, A, which measures the degree of momentum- 
space isotropy of the fluid as follows 

A = ^-l, (2.16) 

PL 

where pr = (T^^ + T2'S')/2 and pi = T^' = -T^ are the effective transverse and longitudinal 
pressures, respectively. If A = 0, the system is locally isotropic. If — 1 < A < the system 
has a local prolate anisotropy in momentum space and if A > the system has a local 
oblate anisotropy in momentum space. In appendix [B] we derive the relation between the A 
parameter defined above and the ^ parameter introduced in Ref. j46[ to quantify the degree 
of local plasma isotropy. For small values of A the relation is A = 4(^/5 + 0{^'^). 

In the 0+1 dimensional model of viscous hydrodynamics one can express the effective 
transverse pressure as pt = p + n/2 and the effective longitudinal pressure as pl = p — U. In 
the case of an ideal equation of state, rewriting fl2.12p in terms of our dimensionless variables 
gives 

At the initial time r = 1, Aq = A(r = 1) is given by 

Ao - I (r^) . (2-18) 



2 VI -SHo/ 

In the limit 11 — 2e/3 we have A —1 and in the limit 11 e/3 we have A oo. 

Positivity of the longitudinal pressure requires A 7^ oo at any time during the evolution 
of the plasma. Note that requiring positivity is a weak constraint on the magnitude of A 
since the formal justification for applying viscous hydrodynamical approximations is the 
neglect of large gradients and higher-order nonlinear terms. This requires that 11 be small 
compared to the pressure, p, i.e. |n| <C p. This can be turned into a quantitative statement 
by requiring that —ap < U < ap, where a is a positive phenomenological constant which 
is less than or equal to 1, i.e. < a < 1. The limit a — > 1 gives the weak constraint of 
—3/4 < A < 00 and for general a requires A_ < A < A_|_ where 

For example, requiring a = 1/3 we would find the constraint — 3/8 < Aq, < 3/4. 



III. APPROXIMATE ANALYTIC SOLUTION OF 0-i-l CONFORMAL HYDRO- 
DYNAMICS 

In this section we present an approximate analytic solution to the 0+1 dimensional con- 
formal 2nd-order hydrodynamical evolution equations. The approximation used will be to 
first exactly integrate the differential equation for the energy density fl2.13ap . thereby ex- 
pressing the energy density as an integral of the shear. We then insert this integral relation 
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into the equation of motion for shear itself (12.13bp and expand in 77. Exphcitly, the solution 
obtained from the first step is 



^-4/3 



1+ / df'ir'V/^Uif' 



(3.1) 



We then solve the second differential equation for 11 approximately by dropping the second 
term in Eq. ( 13. ip and inserting this into the second equation of ( 12.13^ to obtain 



27 cx, 7 k f n + 72 f dril + (72 7 A; f 2 + 96 f^/^) H = 128 . 
This differential equation has a solution of the form 

n ■ ' 



(3.2) 



4/3 



c 




-af2/3j +a(6-l 


) r2/3 


2-6 
V 3 


-af2/3^ 




+ 2G?;° (af2/3 


0,0 y 




-af2/3j 


-G 


(afV3 


b+l\ 
0,1 ^ 



(3.3) 

where iFi is a confluent hypergeometric function, G is the Meijer G function, a = S'jk/ (2c^), 
b = cai^z/ctt, and C is an integration constant which is fixed by the initial condition for 11 at 
f = 1. Requiring n(r = 1) = IIo fixes C to be 



C 



8G?;° (a 


0,0 y 


+ 3 Ho G^;" (^a 


6+1 A 
0,1 j 




[3acA, Ho -8] iFi (^'^ 


— a 


-4a(6-l)iFi 


2-b 
V 3 


— a j 



(3.4) 



To obtain the proper-time evolution of the energy density one must integrate (13.11) using 
(13. 3p . This is possible to do analytically but the answer is rather unwieldy and hence not 
very useful to list explicitly. Below we will use this approximate analytic solution as a 
cross check for our numerics. In the limit f] ^ this solution becomes an increasingly 
better approximation and hence represents the leading correction to ideal hydrodynamical 
evolution in that limit. 

Note that in the limit ca^ — > and c^r — the differential equation above (13.21) reduces 
to an algebraic equation 

nideal Navier-Stokes = 7; 7^ : (3-5) 

97A;r^ 

which, when converted back to dimensionful variables, corresponds to the Navier-Stokes 
solution under the assumption that e = f~^/^. Finally we note that in the large time limit 
Eq. (13.31) simplifies to 



lim n = n 

T— »00 



Ideal Navier-Stokes 



+ e 



-2/3 



(3.6) 



IV. RESULTS 



In this section we present our results of numerical integration of Eq. (I2.13P and present 
consistency checks obtained by comparing these results with the approximate analytic solu- 
tion presented in the previous section. 
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n ' r 



x„ = 0.4 fm/c 



Xq = 1 fm/c 



Xq = 2 fm/c 




FIG. 1: Result for the proper-time evolution of A obtained by numerical integration of Eq. ()2.13p . 
Long-dashed, solid, and short-dashed lines correspond tq = {0.4, 1,2} fm/c, respectively. Transport 
coefficients were the typical strong coupling values given in Eq. (j2.14p . The initial temperature. 
To, is held fixed at Tq = 350 MeV and it is assumed that Aq = for this example. 



A. Time Evolution of A 

Below we present numerical results for the time evolution of the plasma anisotropy pa- 
rameter A. For purpose of illustration we will hold the initial temperature fixed at T = 
350 MeV and vary the starting time tq. This will allow us to probe different values of 
k = TQej/^ = roTo/7 in a transparent manner. Note that, by doing this, each curve corre- 
sponds to a different initial entropy density; however, this is irrelevant for the immediate 
discussion since we are not concerned with phenomenological consequences, only with the 
general mathematical properties of the system of differential equations as one varies the 
fundamental parameters. In Sees. IIV CI and IIVDI we will present the general results as a 
function of the dimensionless parameter k. 



1. Strong Coupling 

In Fig. [1] we show our result for the proper-time evolution of the pressure anisotropy 
parameter. A, obtained by numerical integration of Eq. fl2.13p . The transport coefficients in 
this case are the typical strong coupling values given in Eq. fl2.14p . For purpose of illustration 
we have chosen the initial temperature, Tq, to be held fixed at Tq = 350 MeV and assumed 
that the initial pressure anisotropy, Aq, vanishes, i.e. Aq = 0. 

As can be seen from this figure, when the initial value of the pressure anisotropy is taken 
to be zero it does not remain so. A finite oblate pressure anisotropy is rapidly established 
due to the intrinsic longitudinal expansion of the fluid. Depending on the initial time at 
which the hydrodynamic evolution is initialized, A peaks in the range 0.2 < A < 1. 
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FIG. 2: Comparison of result for A as a function of proper time using numerical integration 
of Eq. (j2.13p and the approximate analytic solution given via Eqs. ()3.3p and (jS.ip . Transport 
coefficients in this case are the typical strong coupling values given in Eq. (j2.14p . The initial 
temperature, To, is taken to be To = 350 MeV, the initial time, to, is taken to be ro = 1 fm/c and 
it is assumed that Ao = for this example. 

2. Comparison with analytic approximation 

As a cross check of our numerical method, in Fig. [2] we compare the result for A obtained 
via direct numerical integration of Eq. f l2.13p and the approximate analytic solution given 
via Eqs. (13.31) and (13. ip . As can be seen from the figure the analytic solution provides a 
reasonable approximation to the true time-evolution of the plasma anisotropy. The param- 
eter A is a particularly sensitive quantity to compare. If one compares the analytic and 
numerical solutions for the energy density, for example, in the strongly-coupled case there 
is at most a 1% deviation between the analytic approximation and our exact numerical in- 
tegration during the entire 10 fm/c of simulation time. Of course, for larger viscosity the 
analytic approximation becomes more suspect but for the weakly-coupled case we find that 
there is at most a 8% deviation between the energy densities obtained using our analytic 
approximation and the exact numerical result. In the limit that fj goes to zero, the analytic 
treatment and our numerical integration agree to arbitrarily better precision. Based on the 
agreement between the two approaches we are confident in our numerical integration of the 
coupled differential equations. 

3. Weak Coupling 

In Fig. [3] we show our result for the proper-time evolution of the pressure anisotropy 
parameter. A, obtained by numerical integration of Eq. (I2.13p . The transport coefficients in 
this case are the typical weak coupling values given in Eq. (I2.15p . For purpose of illustration 
we have chosen the initial temperature. To, to be held fixed at To = 350 MeV and assumed 
that the initial pressure anisotropy, Ao, vanishes, i.e. Aq = 0. 

As can be seen from this figure, as in the strongly coupled case, a finite oblate pressure 
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"1 ' r 



"1 ' r 



— - ^0 = 


0.4 fm/c 


^0 = 


1 fm/c 




2 fm/c 




FIG. 3: Result for the proper-time evolution of A obtained by numerical integration of Eq. ()2.13p . 
Long-dashed, solid, and short-dashed lines correspond tq = {0.4, 1,2} fm/c, respectively. Transport 
coefficients in this case are the typical weak coupling values given in Eq. (j2.15p . The initial 
temperature, Tq, is held fixed at Tq = 350 MeV and it is assumed that Aq = for this example. 



anisotropy is rapidly established due to the intrinsic longitudinal expansion of the fluid. 
In the case of weak coupling transport coefficients a larger pressure anisotropy develops. 
Depending on the initial time at which the hydrodynamic evolution is initialized, A peaks 
in the range 1 < A < 9. 

As can be seen from the tq = 0.4 fm/c result, if the initial simulation time is assumed 
to be small, then very large pressure anisotropies can develop. In that case, in dimensionful 
units, the peak of the A evolution occurs at a time of r ~ 2.3 fm/c. Such large pressure 
anisotropies would cast doubt on the applicability of the 2nd-order conformal viscous hy- 
drodynamical equations, since nonconformal 2nd-order terms and higher-order non-linear 
terms corresponding to 3rd- or higher-order expansions could become important.^ If, in the 
weakly coupled case, the initial simulation time tq is taken to be 0.2 fm/c one would find 
that A would become inflnite during the simulation. This divergence is due to the fact that 
the longitudinal pressure goes to zero and then becomes negative during some period of the 
time evolution. 



B. Negativity of Longitudinal Pressure 

In order to explicitly demonstrate the possibility that A diverges, in Fig.lHwe have plotted 
the evolution the longitudinal pressure over the isotropic pressure {p = e/3), pi/p, obtained 
by numerical integration of Eq. fl2.13p for different assumed initial pressure anisotropies. The 
transport coefficients in this case are the typical weak coupling values given in Eq. (12.151) . 
The initial temperature, Tq, is held fixed at Tq = 350 MeV and it is assumed that tq = 0.2 



^ See Ref. [29j for an example of 2nd-order terms which can appear when conformality is broken. 
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10 12 14 16 18 20 

FIG. 4: Result for the proper-time evolution of the ratio of the longitudinal pressure over the 
pressure, pi/p, obtained by numerical integration of Eq. (j2.13p . Solid, long-dashed, and short- 
dashed lines correspond Aq = {0,-0.5,10}, respectively. Transport coefficients in this case are 
the typical weak coupling values given in Eq. (j2.15p . The initial temperature, Tq, is held fixed 
at To = 350 MeV and it is assumed that tq = 0.2 fm/c for this example. The dotted grey line 
indicates pi = ^ va. order to more easily identify the point in time where the longitudinal pressure 
becomes negative. 



fm/c for this example. 

As this figure shows, if the initial simulation time is too early, the longitudinal pressure of 
the system can become negative. The exact point in time at which it becomes negative de- 
pends on the assumed initial pressure anisotropy. As the initial pressure anisotropy becomes 
more prolate, the time over which the longitudinal pressure remains positive is increased. 
For initially extremely prolate distributions the longitudinal pressure can remain positive 
during the entire simulation time. In the opposite limit of extremely oblate distributions, 
the longitudinal pressure can become negative very rapidly and remain so throughout the 
entire lifetime of the plasma. We note that in the Navier-Stokes limit the initial shear would 
be (no) j^g^^igj. g^Qj^gg = 16r//(9roTo) which, using the initial conditions indicated in Fig. HI 
gives Pl,o/p = —11.1. This means that if one were to use Navier-Stokes initial conditions 
the system would start with an extremely large negative longitudinal pressure. Using tq = 1 
fm/c and Tq =350 MeV improves the situation somewhat; however, even in that case the 
initial Navier-Stokes longitudinal pressure remains negative with Pl,o/p = —1-4. 

What does a negative longitudinal pressure indicate? From a transport theory point 
of view it indicates that something is unphysical about the simulation since in transport 
theory the pressure components are obtained from moments of the momentum-squared over 
the energy, e.g. for the longitudinal pressure 

where /(p) is the one-particle phase-space distribution function. Therefore, in transport 
theory all components of the pressure are positive definite. It is possible to generate negative 
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longitudinal pressure in the case coherent fields as in the case of the early-time evolution 



of the quark-gluon plasma [47|, |48|, |49|, |50|; however, such coherent fields are beyond the 
scope of hydrodynamical simulations which describe the time evolution of a locally color- 
and charge-neutral fiuid. 

This fundamental issue aside, the negativity of the longitudinal pressure indicates that 
the expansion which was used to derive the hydrodynamical equations themselves is break- 
ing down. This expansion implicitly relies on the perturbation described by 11 being small 
compared to the isotropic pressure p. The point at which the longitudinal pressure goes to 
zero is the point at which the perturbation, 11, is equal in magnitude to the background 
around which one is expanding. This means that the perturbation is no longer a small 
correction to the system's evolution and that higher order corrections could become impor- 
tant. Therefore negative longitudinal pressure signals regions of parameter space where one 
cannot trust 2nd-order viscous hydrodynamical solutions. In the following two subsections 
we will make this statement quantitative and extract constraints on the initial conditions 
which allow for 2nd-order viscous hydrodynamical simulation. 



C. Determining the critical line in initial condition space 

For a fixed set of transport coefficients given by {f],c.,t,C\-^} the only remaining freedom 
in the hydrodynamical evolution equations fl2.13p comes from the coefficient 7 (using the 
assumed ideal equation of state) and from the initial conditions through the dimensionless 
coefficient k = Toe^^ and the initial shear IIo. In the next section we will vary these two 
parameters and determine for which values one obtains a solution which, at any point during 
the evolution, has a negative longitudinal pressure. For a given Hq we find that for k below 
a certain value, the system exhibits a negative longitudinal pressure. We will define this 
point in k as the "critical" value of k. Above the critical value of k the longitudinal pressure 
is positive definite at all times. 



1. Strong Coupling 

In Fig.[5]we plot the critical boundary in k (^critical) as a function of the initial value of the 
shear, Uq. Since k is proportional to the assumed initial simulation time Tq increasing k with 
fixed initial energy density corresponds to increasing tq. Assuming fixed initial temperature, 
for an initially prolate distribution, one can start the simulation at earlier times. For an 
initially oblate distribution, one must start the simulation at later times in order to remain 
above the critical value of k. In general, k = tqsI^'^ and our result can be used to set a bound 
on this product. 

In the case of typical strong coupling transport coefficients, the critical value of k at 
Ho = is fccriticai(no = 0) = 0.26. In the case of an ideal QCD equation of state and assuming 
Ho = 0, the constraint is that tq > 7 ^critical TJ,"^, which is numerically tq > 0.14 Tq"^. 
Assuming an initial time of tq = 1 fm/c = 5.07 GeV~^ this implies that Tq > 28 MeV. For 
other initial values of Hq one can use Fig. [5] to determine the constraint. 
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FIG. 5: Critical boundary in k (/^critical) as a function of the initial shear, Ho. Above this line 
solutions have positive longitudinal pressure at all times. Below this line solutions have negative 
longitudinal pressure at some point during the evolution. Transport coefficients in this case are the 
typical strong coupling values given in Eq. (|2.14p . Left limit of plot region corresponds to Aq = — 1 
and right to Aq = oo. 

2. Weak Coupling 

In Fig. [6] we plot the critical boundary in k (/^critical) as a function of the initial value of 
the shear, Hq. Since k is proportional to the assumed initial simulation time tq increasing 
k with fixed initial energy density corresponds to increasing tq. As in the case of strong 
coupling, for an initially prolate distribution, one can start the simulation at earlier times. 
For an initially oblate distribution, one must start the simulation at later times in order to 
remain above the critical value of k. 

In the case of typical weak coupling transport coefficients the critical value of /c at Hq = 
is A;criticai(no = 0) = 0.74. In the case of an ideal QCD equation of state and assuming 
Ho = 0, the constraint is that tq > 7 fccriticaiTj)"^, which is numerically, tq > 0.40 Tq"^. 
Assuming an initial time of tq = 1 fm/c this implies that Tq > 79 MeV. For other initial 
values of IIq one can use Fig. [6] to determine the constraint. 

D. For which initial conditions can one trust 2nd-order viscous hydrodynamical 
evolution? 

As mentioned in Sec. Ill El the requirement that the longitudinal pressure is positive during 
the simulated time only gives a weak constraint in the sense that it merely requires that 
n < p. A stronger constraint can be obtained by requiring instead —ap < 11 < and 
then using this to constrain the possible initial time and energy density which can be used 
in hydrodynamical simulations. In the following subsections we will fix a = 1/3 as our 
definition of what is a "large" correction. For this value of a the initial values of Ho are 
constrained to be between —1/9 < Hq < 1/9. For a given Hq in this range we find that for 
k below a certain value we cannot satisfy the stronger constraint at all simulated times. We 
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FIG. 6: Critical boundary in k (^critical) as a function of the initial shear, Ho. Above this line 
solutions have positive longitudinal pressure at all times. Below this line solutions have negative 
longitudinal pressure at some point during the evolution. Transport coefficients in this case are the 
typical weak coupling values given in Eq. (|2.15|) . Left limit of plot region corresponds to Aq = — 1 
and right to Aq = oo. 

will define this point in k as the "convergence" value of k or /^convergence- Above this value of 
k = /^convergence the sheai Satisfies the constraint —p/3 < 11 < p/3 at all simulated times and 
therefore represents a "reasonable" simulation. 

1. Strong Coupling 

In Fig. [7] we plot the " COIlVCrgGIlCG boUIlcicLry in k (^convergence) ^ 

function of the initial 

shear, IIq. In the case of typical strong coupling transport coefficients the convergence value 
of /c at Ho = is /cconvergence(no = 0) = 1.58. In the case of an ideal QCD equation of state 
and assuming IIq = 0, the constraint is that tq > 7 /^convergence ^o~^ which is numerically 
To > 0.85 Tq^. Assuming an initial time of tq = 1 fm/c this implies that Tq > 167 MeV. For 
other initial values of Uq one can use Fig. [7] to determine the constraint. 

2. Weak Coupling 

In Fig. [8] we plot the "convergence boundary" in k (/Cconvergence) as a function of the initial 
shear. Ho. In the case of typical weak coupling transport coefficients the convergence value 
of /;; at Ho = is /i^convergencelHo = 0) = 10.9. In the case of an ideal QCD equation of state 
and assuming Ho = 0, the constraint is that tq > 7 /^convergence TJ^"^, which is numerically 
Tq > 5.9 Tq"^. Assuming an initial time of tq = 1 fm/c = 5.07 GeV~^ this implies that 
To > 1.16 GeV. For other initial values of IIq one can use Fig. [H]to determine the constraint. 
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FIG. 7: Convergence boundary in k (^convergence) as a function of the initial shear, XIq. Above this 
hue solutions satisfy the convergence constraint. Transport coefficients in this case are the typical 
strong coupling values given in Eq. (|2.14p . 

E. What does this imply for higher dimensional hydrodynamical simulations? 

If one proceeds to more realistic simulations in higher dimensional boost invariant treat- 
ments, e.g. 1+1 and 2+1, the spatial variation of the initial conditions and time evolution 
in the transverse plane have to be taken into account. In addition, new freedoms such 
as the initial fluid flow field and additional transport coefficients arise; however, to first 
approximation one can treat these higher dimensional systems as a collection of 0+1 dimen- 
sional systems with different initial conditions at each point in the transverse plane. Within 
this approximation one would quickly find that there are problems with the hydrodynamic 
treatment at the transverse edges of the simulated region. 

This happens because as one goes away from the center of the hot and dense matter, 
the energy density (temperature) drops and, assuming a fixed initial simulation time tq, one 
would find that at a finite distance from the center the condition k > ^critical would be violated 
by the initial conditions. In these regions of space, hydrodynamics would then predict an 
infinitely large anisotropy parameter. A, casting doubt on the reliability of the hydrodynamic 
assumptions. Even worse is that at a smaller distance from the center one would cross 
the "convergence boundary" in k, /^convergence , and therefore not fully trust the analytic 
approximations used in deriving the hydrodynamic equations (conformality, truncation at 
2nd order, etc.). 

Of course, an approximation by uncoupled 0+1 systems with different initial conditions 
would not generate any radial or elliptic flow; however, we flnd empirically that the picture 
above holds true in higher-dimensional simulations, justifying the basic logic. For exam- 
ple, using strongly-coupled transport coefficients and assuming an initially isotropic plasma 
(Ho = 0), we found in Sec. IIV C II that ^critical = 0.26. In terms of the initial temperature 
this predicts that when starting a simulation with tq = 1 fm/c, one will generate negative 
longitudinal pressures for any initial temperature Tq < 28 MeV. 

We will now compare this prediction with results for the longitudinal pressure extracted 
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FIG. 8: Convergence boundary in k (Aiconvergence) as a function of the initial shear, XIq. Above this 
hue solutions satisfy the convergence constraint. Transport coefficients in this case are the typical 
weak coupling values given in Eq. (|2.15p . 

from the 2+1 dimensional code of Luzum and Romatschke gim. In Fig.i we show fixed 
r snapshots of the longitudinal pressure. The runs shown in Fig. [9] were performed on a 
69^ transverse lattice with a lattice spacing of 2 GeV^^ using Glauber initial conditions 
starting at tq=1 fm/c, an initial central temperature of Tq = 350 MeV, zero initial shear 
and zero impact parameter. For these runs we have used the realistic QCD equation of 
state used in Ref . . In the left panel of Fig. [H] the transport coefficients were set to the 
typical strong coupling values given in Eq. (12.141) . except with = due to the fact 
that the code used did not include this term in the hydrodynamic equations. Based on the 
initial transverse temperature profile and our estimated critical initial temperature, in the 
strong-coupling case we expect negative longitudinal pressures to be generated at transverse 
radius r > 10 fm. As can be seen from the left panel of Fig. [9l at the edge of the simulated 
region the longitudinal pressure becomes negative starting already at very early times. The 
transverse radii at which this occurs is in good agreement with our estimate based on the 
0+1 dimensional critical value detailed above. 

Based on our convergence criterium detailed in Sec. lIVDI we found, in the strong-coupling 
case, that /cconvergence(no = 0) = 1.58. Assumiug To = 1 fm/c this translates into a minimum 
initial temperature of 167 MeV. Based on the transverse temperature profile used in the 
run shown in the left panel of Fig. [9] this results in a maximum transverse radius r ~ 6.8 
fm. At radii larger than this value it is possible that higher order corrections are large 
and therefore the applicability of 2nd-order viscous hydrodynamics becomes questionable. 
Since this temperature is greater than the typical freeze-out temperature used, Tf ~ 150 
MeV, this means that in the strong coupling limit it is relatively safe to use hydrodynamical 
simulations. However, one should be extremely careful with the transverse edges. 

The situation, however, is not as promising in the weak-coupling case. To see this ex- 
plicitly, in the right panel of Fig. [9] we show the longitudinal pressure resulting from a run 
with weak coupling transport coefficients (12.150 . Based on the initial transverse temperature 
profile and our estimated critical initial temperature, in the weak-coupling case we expect 
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FIG. 9: Evolution of tlie longitudinal pressure in proper-time obtained from the 2+1 dimensional 
viscous hydrodynamics code of Ref. [5]. Horizontal axis is the distance from the center of the 
simulated region. In the left panel we show the result obtained using the typical strong coupling 
values given in Eq. (j2.14p but with c^^ = 0. In the right panel we show the result obtained using the 
typical weak coupling values given in Eq. (j2.15p but with cx^ = 0. The runs shown used Glauber 
initial conditions with an initial central temperature of Tq = 350 MeV, initial time tq = 1 fm/c 
and n;;(ro) = 0. 



negative longitudinal pressures to be generated at transverse radius r > 8 fm. Comparing 
this prediction to the results shown in the right panel of Fig. [9] we see that the situation is 
even worse than expected. By the final time of 4.5 fm/c the entire central region has very 
low or negative longitudinal pressure. We note that at that time the radius at which the 
temperature has dropped below the freeze-out temperature is around 7.3 fm so the region 
where the longitudinal pressure is negative (or almost negative) is still in the QGP phase. 

In terms of convergence, we remind the reader that based on our convergence criterium 
detailed in Sec. IIVDI we found that in the weakly-coupled case /cconvergence(no = 0) = 10.9. 
Assuming tq = 1 fm/c we found that the initial central temperature should be greater 
than 1.16 GeV. As can be seen in Fig. [H] the corrections to ideal hydrodynamics are sizable 
so this again points to the possibility that there are large corrections to the 2nd-order 
hydrodynamic equations. Based on this, it would be questionable to ever apply 2nd-order 
viscous hydrodynamics to a weakly-coupled quark-gluon plasma generated in relativistic 
heavy-ion collisions. At the very least one would need to include nonconformal 2nd-order 
terms and 3rd-order terms in order to assess their impact. 



V. CONCLUSIONS AND OUTLOOK 

In this paper we have derived two general criteria that can be used to assess the applica- 
bility of 2nd-order conformal viscous hydrodynamics to relativistic heavy-ion collisions. We 
did this by simplifying to a 0+1 dimensional system undergoing boost invariant expansion 
and then (a) requiring the longitudinal pressure to be positive during the simulated time or 
(b) requiring a convergence criterium that |n| < p/3 during the simulated time. We showed 
that these requirements lead to a non-trivial relation between the possible initial simulation 
time To, the initial energy density eo, and the initial value of the fiuid shear tensor. Ho. As 
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a cross check of our numerics we presented an approximate analytic solution of 2nd-order 
conformal viscous hydrodynamical evolution which represents the leading correction to 0+1 
dimensional boost-invariant ideal hydrodynamics in the limit ?7/s — > 0. 

The constraints derived here were then shown to provide guidance for where one might 
expect 2nd-order viscous hydrodynamics to be a good approximation in higher-dimensional 
cases. We found that the prediction of our criticality bound was in reasonable agreement 
with where the longitudinal pressure becomes negative in 2+1 dimensional viscous hydro- 
dynamical simulations. Based on these findings it seems possible to estimate where one 
obtains convergent/trustable 2nd-order viscous hydrodynamical simulations based solely on 
the initial conditions and analysis of the hydrodynamical evolution equations themselves. 

In closing we mention that another outcome of this work is that we have shown that 
it is possible to use hydrodynamical simulations to predict the proper-time dependence 
of the plasma momentum-space anisotropy as quantified by the A or ^ parameters. This 
can be used as input to calculations of production of electromagnetic radiation from an 
anisotropic plasma [9|,0,|h, T^, 13, 14, 1^, calculations of quarkonium binding/polarization 



in anisotropic plasma [52, l53| , and also to assess the phenomenological growth rate of plasma 
instabilities on top of the mean colorless fluid background (see Ref. [s^] and references 
therein). The findings here present a complication in this regard since phenomenological 
studies will require knowledge of A in the full transverse plane. As we have shown, 2nd- 
order hydrodynamical simulations predict that this parameter can become infinite in certain 
regions. In these regions one would no longer trust the predictions of the hydrodynamical 
model and additional input would be required. 
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APPENDIX A: NOTATION AND CONVENTIONS 

We summarize the conventions and notation we use in the main body of the text: 

• The metric for a Minkowski space in the curvilinear coordinates {T,x,y,() is Qfiu 

di&g{grr, Qxx, 9yy, 9cc) = (1> "l, "l, -^^) • 

• A'^'^ = g^'^ — u^u^ is a projector orthogonal to the fluid velocity, u^A'^'^ = 0. 

• The comoving time derivative: D = u^Dc,. 

• The comoving space derivative: = A^'^D^. 
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• The brackets ( ) denote an operator that is symmetric, traceless, and orthogonal to 
the fluid velocity: 

A^,B,) = (a;A^ + A^A^ - ^ A"^A,.^ A^B^ . (Al) 

• The symmetric and anti-symmetric operators: 

= ^ {A,,B, + A,B^) , (A2) 
= ^ {A^B, - A,B^) . (A3) 

APPENDIX B: RELATION BETWEEN A AND ^ 

In this appendix we derive the relation between the anisotropy parameter A introduced 
in this paper and the ^ parameter introduced in Ref. [46|. In the general case ^ is defined by 
taking an arbitrary isotropic distribution function /iso(p) and stretching or squeezing it along 
one direction in momentum space to obtain an anisotropic distribution. Mathematically this 
is done by introducing a unit vector n which defines the direction of anisotropy, an anisotropy 
parameter — 1 < ^ < cxd, and requiring /(p) = /iso^\/p^ + ^(p " fi)^j • Fixing fi = z to define 
the longitudinal direction and assuming massless particles, is it straightforward to evaluate 
the transverse and longitudinal pressures through the components of the stress-energy tensor 

PT = I (T-- + T-) = lj-^s ( v'i^TM) , (Bl) 

and 

By a change of variables to p = a/p^ + ^pl and the use of spherical coordinates one can 
show that 



4f V ' VI J 

and 



PL 



where and are the isotropic transverse and longitudinal pressures which are obtained 
from /iso, respectively. Combining the above relations and using = = e^'^°/3, where 
e'^° is the isotropic energy density, we obtain the following expression for A 

A = i«-3)+f((l + 0^-l)". (B6) 

In the small ^ limit 

hmA = ^C + 0{e), (B6) 
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and in the large ^ limit 

\imA = h + 0{^). (B7) 
For general ^ one needs to invert flB5l) numerically in order to obtain ^ as a function of A. 
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